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We present the results of extensive Monte Carlo simulations of Ising models with algebraically 
decaying ferromagnetic interactions in the regime where classical critical behavior is expected for 
these systems. We corroborate the values for the exponents predicted by renormalization theory 
for systems in one, two, and three dimensions and accurately observe the predicted logarithmic 
corrections at the upper critical dimension. We give both theoretical and numerical evidence that 
above the upper critical dimension the decay of the critical spin-spin correlation function in finite 
systems consists of two different regimes. For one-dimensional systems our estimates for the critical 
couplings are more than two orders of magnitude more accurate than existing estimates. In two and 
three dimensions we give, to our knowledge, the first results for the critical couplings. 
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I. INTRODUCTION 

The critical behavior of Ising models with long-range interactions has attracted much Attention during the last 
three decades. For the one-dimensional case, some analytical results have been obtakiedjEfllJ as well as a number 
of numerical results. jXhc-pumerical results apply to both inverse-square interactionstjcj and geaecal algebraically 
decaying interactions .c3~c3 Special mention deserves the work by Anderson, Yuval, and HamannjHa~E!l which greatly 
* ■ stimulated the interest in spin chains with long-range-interactions. They also developed a renormalization-like ap- 
proach to the one-dimensional inverse-square model.E2rEj Further renormalization-group studies of this_particular case 
are presented in Refs. |T^]3^ -|34|. A major contribution was made by Fisher, Ma, and NickeO and Sak,E3 who obtained 
' renormalization predictions for the critical exponents of models-pf general dimensionality d < 4 with algebraically 
decaying interactions (obtained independently by Suzuki et aZEH). Other works concerning d > 1 are two conjec- 
tures on, respectively, the boundary between long-range and short-range behavior and the boundary between classical 



(mean-field) and nonclassical behavior, both by Stelljla a (refuted) conjecture by Grimthsjia a rigorous confirmation 
i of the upper critical dimension by Aizenman and Fernandez and a variational approach to the Ising model with 
t-H ' long-range interactions.!^ Furthermore, Monte Carlo Emulations have been carried out for one particular choice of 
the spin-spin interaction in a two-dimensional model.E-U However, to our knowledge, neither any further verifications 
of the renormalization predictions nor any other results are available for higher-dimensional (d > 1) models. To 
conclude this summary, we-mention thatp-the. one-dimensional g-state Potts model with long-range interactions has 
been studied analytically,BO numericallyJl3'EJ and in a mean-field approximation on the Bethe lattice.cil 
— «, | Why are these models interesting? In the first place from a fundamental point of view: They enable us to study the 
■ influence of the interaction range on the critical behavior. E.g., in one-dimensional systems long-range order is only 
possible in the presence of spin-spin interactions which decay sufficiently slowly. In the borderline (inverse-square) 
I . case, the ID model displays a remarkable behavior: At the critical temperature the order parameter exhibits a finite 
jump (see Sec. ||), but the free energy has an essential singularity such that all thermal properties are smooth. Ir 
sense, the phase transition can be regarded as the one-dimensional analog of a Kosterlitz-Thouless transition! 
although the jump in the magnetization is not present there, as follows from the Mermin- Wagner theoremJfD Just 
as d = 2 is the lower critical dimension for the two-dimensional XY model with short-range interactions, a = 1 is 
a critical decay rate in a one-dimensional system with interactions decaying as r~( 1+CT ), see Ref. [32[ With respect 
to higher-dimensional systems, we note that the decay rate of van der Waals forces in realistic three-dimensional 
systems is only slightly faster than at the boundary between short-range (Ising-like) and long-range critical behavior. 
The question of criticality in ionic systems, where the (screened)— Coulomb interactions might lead to effectively 
algebraically decaying interactions, appears still open to debate.E3~E5l It has also been claimed that-exponents in the 
long-range universality class have been observed experimentally in a ferromagnetic phase transitionEI Recently, it has 
been derived that critical fluctuations may give rise to long-range Casimir forces (decaying much more slowly than 
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van der Waals interactions) between uncharged particles immersed in a critical fluid.Ea Furthermore, it was shown by 
Anderson and YuvalE3£3 that the Kondo problem corresponds to a one-dimensional Ising model with a combination 
of inverse-square and nearest-neighbor interactions. Yet another application follows from Ref. 22, where it was shown 
that random exchange (Levy-flight) processes can generate effective interactions which decay algebraically. Hence, the 
universal critical properties of the nonequilibrium steady state of these systems are those of the long-range equilibrium 
Ising models studied in this paper. Finally, the realization that the upper critical dimension can be varied by tuning 
the decay rate of the interaction led to a special application of these models in Ref. |5^. Here, they were used to 
analyze a long-standing controversy on the universality of the renormalized coupling constant above the upper critical 
dimension. 

In this article, we present accurate numerical results for Ising systems with algebraically decaying interactions 
in one, two, and three dimensions. Until now, the long-range character of the spin-spin interactions has been the 
main bottleneck for the examination of these systems by means of numerical methods (and, in fact, also for their 
analytical solution). All previously published numerical results therefore rely cp, various extrapolations based on 
data for small systems. However, the advent of a novel Monte Carlo algorithmEJ for the first time enabled us to 
efficiently simulate these systems. The high accuracy of the results opens several perspectives: i) verification of the 
renormalization predictions for the critical exponents; ii) accurate observation of logarithmic corrections at the upper 
critical dimension; iii) first estimates of the critical temperatures of two- and three-dimensional systems with long- 
range interactions; iv) verification of previously obtained estimates of the critical temperatures of one-dimensional 
systems, which in addition implies a check on the various extrapolation methods that have been developed; v) 
verification of predicted bounds on the critical temperatures; vi) verification of a conjecture on the behavior of the 
critical temperature as a function of the decay parameter. Another problem one encounters in the simulations is 
the large parameter space: The simulations for a set of different temperatures and system sizes have to be repeated 
for a range of values of the decay parameter and for d = 1, 2, 3. The total computing time dedicated to the results 
presented in this paper amounts to approximately two CPU-years on a modern workstation. 

The outline of this paper is as follows. In Sec. 0, we sum up the known rigorous results for the Ising chain with 
long-range interactions. In Sec. HI, we review the renormalization scenario of these models and derive the finite-size 
scaling behavior of several quantities. These include the corrections to scaling, both at and above the upper critical 
dimension. Our numerical results are presented and analyzed in Sec. |lV] and compared with previously obtained 
results. Finally, we summarize our conclusions in Sec. The Appendix contains technical details concerning the 
application of the long-range Monte Carlo algorithm to the models studied in this paper. 



II. RIGOROUS RESULTS FOR THE ONE-DIMENSIONAL CASE 



For the one-dimensional case, the Hamiltonian is given by 

H = Y^J(i-j)s i 8 j , (1) 

y 

where the sum runs over all spin pairs. We are particularly interested in algebraically decaying interactions, i.e, 
J(n) oc n~ a . To ensure that the energy of the system does not diverge, it is required that a > 1. In 1968, RuelleEl 
rigorously proved the absence of long-range order in a spin chain with ferromagnetic spin-spin couplings J(i — j) such 
that the sum 

N 

71=1 

does not diverge in the limit N — > oo. For. algebraically decaying interactions, this implies the absence of a phase 
transition for a > 2. Shortly later, DysonB proved the existence of a phase transition if the sums , J(n) and 
^^ =1 (loglogn) [n 3 J(n)] 1 both converge, for positive and monotonically decreasing J(n). In particular, a phase 
transition occurs for J(n) oc n~ a with I < a < 2. This partly corroborated the conjecture of Kac and Thompson,E3 
viz. that there is a phase transition for I < a < 2. Furthermore, DysonO was (as were — much later — also Rogers 
and Thompsoncl) able to replace Ruelle's condition with a stronger one, which however still left the case a = 2 
undecided. This alaa holds for an even more stringent criterion by Thoulessp who generalized the argument of 
Landau and LifshitzE3 for the absence of a phase transition in an Ising chain with short-range interactions. However, 
Thouless argued on entropic grounds that if a phase transition exists for a = 2, the magnetization must have a 
discontinuity at the transition point. This was later dubbed the "Thouless effect" by Dyson, who proved it to occur 
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in the closely related hierarchical model. Simon and Sokal made Thouless' argument partially rigorous,Q but later 
Aizenman et a/El showed that, although a discontinuity in the order parameter is indeed present if there is a phase 
transition, his argument does not account for this. Namely, Thouless had assumed that the spin-spin correlation 
function (soS r ) — (so)(s r ) vanishes in the limit r — > oo, whereas actually the critical exponent r\ is equal to 1 in this 
case. Meanwhile, Frohlich and SpencerEI had been able to rigorously prove the existence of a phase transition in the 
borderline case and thus to corroborate the Kac-Thompson conjecture for a — 2 as well. Another interesting point 
is the rigorous proof for the existence of an intermediate ordered phase in the one-dimensional model with inverse- 
square interactions, where the two-point correlation function exhibits powetrlaw decay with an exponent which varies 
continuously in a finite temperature range below the critical temperature.llil 



III. FINITE-SIZE ANALYSIS OF THE CRITICAL BEHAVIOR 



Already in a very early stage of the history of the e-expansion, Fisher, Ma, and Nickel analyzed the critical behavior 
of d-dimensional systems (d < 4) with long-range interactions decaying as r~( d+a \ with a >-|0.E3 They concluded 
that the upper critical dimension is ffiiyen by d u — 2a, as was previously conjectured by StellE3 and later rigorously 
proven by Aizenman and Fernandez. El For more slowly decaying interactions, < a < d/2, the critical behavior is 
classical, whereas the critical exponents assume_nonclassical, continuously varying values for d/2 < a < 2. For a > 2 
they take their short-range (Ising) values. Sak,E3 however, found that already for a > 2 — rj sr the critical behavior is 
Isinglike, where 77 sr denotes the exponent 77 in the corresponding model with short-range interactions. In this article 
we concentrate on the classical range, for which we have performed extensive. Monte Carlo simulations of spin models 
in d = 1, 2, 3. The nonclassical range will be the subject of a future article.Ea 

We briefly outline the renormalization scenario for these models, in order to derive the finite-size scaling rela- 
tions required to analyze the numerical data. We start from the following Landau-Ginzburg- Wilson Hamiltonian in 
momentum space, 

H((j)k)/kBT= i^(j CT fc CT +j 2 fc 2 +ro)0k0-k + ^X!X!X!^ k i^ k 2 < ? !,k 3^~ki-k2~k3 ~hJ~^<t>k=o ■ (3) 

k ki k2 k3 

The jo-fc' 7 term arises from the Fourier transform of the interactions decaying as r^( d+a \ The j2k 2 term normally 
representing the short-range interactions-is included because it will appear anyway in the renormalization process and 
will compete with the long-range term.Lj Under a renormalization transformation with a rescaling factor b = e l , the 
term j a k a is transformed into j a k /(T , with k' = k&. To keep the coefficient of the k a term fixed, we rescale the field 
0k to <j)' k , — b~ rT / 2 (f>\ s _. Thus, the coefficient of the k 2 term decreases as b"~ 2 and the coefficient of the </> 4 term changes 
proportional to 6 2<T ~ d . Hence, the Gaussian fixed point dominates the renormalization flow for a < d/2, which is the 
situation studied in this paper. 

For the sake of generality we treat here the case of an n-component order parameter with O(n) symmetry. The 
renormalization equations are then given by 

= ar a + a ( n + 2 M C - r o) i (4a) 

^ = eu - a(n + 8)w 2 , (4b) 
dl 

where (n + 2) and (n + 8) are the usual factors arising from the tensorial structure of the interaction part of the 
Hamiltonian and e = 2a — d. These equations are not complete to second order, because the 0(u 2 ) term is missing 
in Eq. @. 

We first consider the case e < 0. The solution of the second equation is given by 

u{l) = ue £l — 1 , (5) 

1 + S £(£±8)( e e/_1) 

where u denotes the value of u at I = 0. This yields, to leading order in u, the following solution for the first equation, 

(n+2)/(n+8) 



ro{l) = [ro + ac(n + 2)u/(d - a)} e al 



1 



1 + ^+^u{e £l - 1) 



ac{n + 2)ue ei /(d - a) 

, . o(n+8)-/ s r 7T ' W 

1 H — L - — -u(e E1 — 1) 
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with fo = i'o(l = 0). The first factor between square brackets is proportional to the reduced temperature t = 
(T — T c )/T c and the last term is the so-called shift of the critical temperature. The factors [1 + a(n + 8)u(e el — 
m Eqs. (§) and (|) are higher-order corrections in u. Under successive renormalization transformations, u approaches 
the value u* — and the Gaussian fixed point (0,0) is thus indeed stable. The pertinent renormalization exponents 
are: yt = a, = (d + cr)/2, and j/i = 2a — d. 



At e = 0, the Gaussian fixed point becomes marginally stable. Solving Eq. (4b) leads to 



1 + a(n + 8)ul 



(7) 



where the superscript "uc" indicates that we are operating at the upper critical dimension. This solution can be used 
to solve, again to leading order in u, Eq. (fla|), yielding 



*o c (0 = [fo + ac{n + 2)ft/(d/2)] e 
or, in terms of the rescaling factor b, 

rg" = [f + ac(n + 2)u/(d/2)]b a 



1 



1 + a(n + 8)ul 



(n+2)/(n+8) _ ac(n + 2)u/(d/2) 
1 + a(n + 8)ul 



1 



1 + a(n + 8)12 In 6 



(ra+2)/(n+8) 



ac(n + 2)u/{d/2) 
1 + a(n + 8)12 In 6 



(8) 



(9) 



Since a is fixed at d/2 the factor d/2 in the last term is identical to the corresponding factor (d — a) in Eq. (||). 
Further comparison of Eqs. (||) and (||) shows that above the upper critical dimension the leading shift of the critical 
temperature is proportional to 6 e , whereas this factor vanishes at the upper critical dimension itself and the factor 
{e el — l)/e in the second-order correction turns into a ln& term, yielding a logarithmic shift of the form 1/ (A\n.b + B). 

From the solutions of the renormalization equations we can derive the scaling behavior of the free energy and of 
(combinations of) its derivatives. For the case e < the free energy density / scales, to leading order, as 



f(t,h,u,l/L) = b- d f (b yt [t + aub yi - yt ] ,b Vi >h,b v >u,b/L) +g 



(10) 



where a — —ac(n + 2)/(d — a) and we have included a finite-size field L . g denotes the analytic part of the 
transformation. We abbreviate the first term on the right-hand side as b~ d f(t' , hi ,u' ,b/ L). However, we must take 
into account the fact that, for T < T c , the free energy is singular at u = 0. This makes u a so-called dangerous 
irrelevant variable; see, e.g., Ref. |5^. As discussed in Ref. |5^, the correct finite-size scaling properties are obtained by 
setting b = L and making the substitution <j)' = 0/ii' 1 / 4 . This leads to a new universal function, /, with 



f(t',h',u',l)+g = f(t,h) , 



(11) 



where t = t' /u' 1 / 2 and h = h'/u' 1 / 4 . The analytic part of the transformation also contributes to the singular 
dependence of the free energy on t (see, e.g., Ref. |6(], Ch. VI, § 3): Despite the regularity of this term in each single 
renormalization step, the infinite number of steps still leads to the build-up of a singularity. This contribution, denoted 
by g, is absorbed in / as well. Setting b = L and combining Eqs. (|l0|) and (JTTI) yields 



/ t,h,u, 



= L- d f L 



,1/2 



[t + auL Vi - '' 



J /ill. 



/4 



(12a) 
(12b) 



Here, we have introduced the exponents y* = yt — yi/2 — d/2 and = j/h — y-jA = id/ A. The corresponding critical 
exponents indeed assume their fixed, classical values; a = 0, /3 = 1/2, 7 = 1, 5 = 3. The exponent 7 is singled out 
here as a special case; even without taking into account the modification of yt and ?/h due to the dangerous irrelevant 
variable one obtains the classical value 7=1. Since the correlation length exponent v = l/y t (it is not affected by 
the singular dependence of the free energy oruu), we see that hyperscaling is violated, which is a well-known result for 
systems above their upper critical dimension.^ The rescaling of the pair-correlation function g (decaying proportional 
to \/r d ~ 2+ri ) relates the exponent rj to the rescaling factor of the field, yielding r\ = 2 — a. Note that this contrasts 
with the short-range case (a = 2), where 77 assumes its mean- field value for all dimensionalities d > 4. This implies 
that direct experimental measurement of either v or 77 offers a way to discern whether the interactions in a system are 
mean-field-like (y = 1/2, r/ — 0) or have the form of a slowly decaying power-law. Below the upper critical dimension, 
however, the finite-size scaling behavior of the spin-spin correlation function is (apart from a volume factor) identical 
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to that of the magnetic susceptibility x- This relation yields a contradiction above the upper critical dimension, 
since x depends on the scaled combination tL Vt , instead of tL Vt . Indeed, the susceptibility diverges as i~ 7 and the 
finite-size behavior of x is thus xl oc L 1Vt = L d f 2 , corresponding to g oc L~ d / 2 . On the other hand, if one assumes 
that the finite-size behavior of the correlation function is identical to its large-distance behavior, one expects that 
g oc L~( d ~ 2 + ? ?) = L~( d ~ a \ Only at the upper critical dimension, d u = 2er, these two predictions coincide. We will 
return to this point at the end of this section. Furthermore, we will examine the behavior of the spin-spin correlation 



function in Sec. IV 



At the upper critical dimension itself, i.e. at e = 0, the free energy density scales as 



/ ( *, h, u, — 



= b- d f 



L- d f 



(l+/3 M ln6)(™+ 2 )/("+ 8 ) 



t + ab~ yt 
1 



(1 + / 3 u lni)(»+2)/(«+8)-i/2 «V3 



'1 + (3u\nbf/( n +%) 



,b Vk h, 



1 + j3u In b L 



T I + // 



(l+/3ulnL) 6 /("+ 8 ) 



,^ h ^7j[l + /^mL] 1/4 



(13a) 
(13b) 



where (3 = a(n + 8) and we have set b = L in the last line, u is now a marginal variable and although we again have 
to perform the substitution <fi — > <f>' (the Gaussian fixed point is marginally sta ble) , the e xpon ents yt and y^ coincide 
with yl and y^, respectively, because y\ vanishes. Thus, the scaling relations (12b) and ( |l3b| ) differ to leading order 
only in the logarithmic factors arising in the arguments of /. 

As usual, the finite-size scaling relations are now found by taking derivatives of the free energy density with respect 
to the appropriate scaling fields. In the Monte Carlo simulations we have sampled the second and the fourth moment 
of the magnetization density, the dimensionless amplitude ratio Q = (m 2 ) 2 / (m A ) (which is directly related to the 
Binder cumulantEll) , and the spin-spin correlation function over half the system size (for even system sizes). The 
second moment of the magnetization density is (apart from a volume factor) equal to the second derivative of the free 
energy density with respect to h, 



(m 2 ) = L 



- d d 2 f 



2rf, 



-1/2^(2) I L y' t 



,1/2 



'1/4 



(14) 



where stands for the second derivative of / with respect to its second argument and i = t + auL Vi yt . At e = 0, 
logarithmic factors do arise not only in the arguments of f^ 2 \ but also in the prefactor, 



{m2) = L2yb ^ d fl + /3ulnL 



1/2 



(2) 



(1 + pu\iiLY n + 2 )/( n + s )- 1 / 2 u 1 / 2 



t + aL~ 



(1 + /3ulnL) 6 /("+ 8 ) 



^^-^[l + pulnL] 1 /*) . (15) 



For the fourth magnetization moment similar expressions hold and in the amplitude ratio Q all prefactors divide out, 
both for e < and e = 0. Thus we find that the ratio Q is given by a universal function Q, 



Ql(T) = Q[L Vi — 



t 



1/2 



(16) 



where we have omitted the h dependence of Q, since we are only interested in the case_i = 0. The additional term 
proportional to qi arises from the h dependence of the analytic part of the free energyo and the ellipsis stands for 
higher powers of L d ~ 2Vh (faster-decaying terms). At e = 0, t must be replaced by the first argument within square 
brackets in Eq. ( |l3b|) , multiplied by the factor (1 + /3ulni) 1 / 2 -(™+ 2 )/("+ 8 ). Finally, we may derive the finite-size 
scaling behavior of the spin-spin correlation function g(r) by differentiating the free energy density to two local 
magnetic fields, which couple to the spins at positions and r, respectively, and assuming that the finite-size behavior 
is identical to the r dependence of g. If we do not take into account the dangerous irrelevant variable mechanism, we 
find g oc L 2yh ~ 2d = L~( d ~ <J ), just as we found before from r/ = 2 — a. However, replacing j/h by t/j* yields g oc L~ d / 2 , 
in agreement with the L dependence of the magnetic susceptibility. This clarifies the difference between the two 
predictions: At short distances (large wave vectors), the j a k a (f>k4>-k term will be the dominant term in the Landau- 
Ginzburg- Wilson Hamiltonian and there is no "dangerous" dependence on u. Hence, the finite-size behavior of the 
spin-spin correlation function will be given by L~( d ~ 2+, 7). For k = 0, the coefficient of the <p 2 term vanishes and thus 
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the ucj) 4 term is required to act as a bound on the magnetization. To account for this singular dependence on u, we 
rescale the field, which implies that y^ is replaced by y^ and g scales as L 2Vh ~ 2d . In a finite system, the wave vectors 
assume discrete values, k = (n x ,n y , n z )2ir/L, and thus it is easily seen that even for the lowest nonzero wave vectors 
jo-fc^^k^-k constitutes the dominant bounding term on the magnetization. Namely, the coefficient of the <j> term 
contains a volume factor L~ d [cf. Eq. &\ and this term is thus (above the upper critical dimension) a higher-order 
contribution decaying as L 2<J ~ d . 

IV. NUMERICAL RESULTS AND COMPARISON WITH EARLIER RESULTS 

A. Simulations 

We have carried out Monte Carlo simulations for systems described by the Hamiltonian 

n/k B T = -J2j(,h-r j \)s i s j , (17) 

«> 

where the sum runs over all spin pairs and periodic boundaries were employed. The precise form of the (long-range) 
spin-spin interaction J(r) as used in the simulations was chosen dependent on the dimensionality. For d — 1 we have 
followed the conventional choice J(r) — K/r d+a (with discrete values for r), as this allows us to compare all our 
results (including nonuniversal quantities) to previous estimates. However, as explained in Ref. |54] and the Appendix, 
this discrete form requires the construction of a look-up table, which becomes inefficient for higher dimensionalities. 
For d = 2 we have thus applied an interaction which is the integral of a continuously decaying function, 

J(|r|) = K f X + 2 dx f V+2 dy r -^ d+ ^ , (18) 

where r = (r x ,r y ) and r = |r|. In d = 3 the corresponding volume integral was used for J(r). This modification 
of the interaction does only change nonuniversal quantities like the critical temperature, but should not influence 
the universal critical properties like the critical exponents and dimensionless amplitude ratios, since the difference 
between the continuous and the discrete interaction consists of faster decaying terms that are irrelevant according to 
renormalization theory. Details concerning the simulations can be found in the Appendix. 

The following system sizes have been examined: chains of length 10 < L < 150000, square systems of linear size 
4 < L < 240, and cubic systems of linear size 4 < L < 64. At the upper critical dimension simulations for even larger 
systems have been carried out in order to obtain accurate results from the analyses: L = 300000 in d = 1 and L = 400 
in d = 2. (I.e., in terms of numbers of particles the largest system size for d = 2 is considerabrw-smaller than for 
d = 1 and d = 3.) For the simulations we used a new cluster algorithm for long-range interactions. til This algorithm 
is 0{L d+z ) times faster than a conventional Metropolis algorithm, where z is the dynamical critical exponent. For 
systems displaying mean-field-like critical behavior, we expect z = d/2 and the efficiency gain in our simulations is 
thus of the order of 10 8 for the largest system sizes. For each data point we have generated between 10 6 and 4 x 10 6 
Wolff clusters. 



B. Determination of the critical temperatures, the amplitude ratio Q, and the thermal exponent 

The critical couplings K c of these systems have been determined using an analysis of the amplitude ratio Q. The 
finite-size scaling analysis was based on the Taylor expansion of Eq. ([if]), which for e < reads: 

Ql(T) = Q+ Pl iL v ' + p 2 t 2 L 2 y' + P3 PL 3yt +■■■ + qi L d - 2 ^ + ■■■ + q 3 L Vi + ■■■ . (19) 

The term proportional to a in i yields a contribution q2L Vi ^ 2 = q2L a ~ d l 2 and the term q 3 L Vi comes from the 
denominator in Eq. (^). The coefficients pi and qi are nonuniversal. In addition to the corrections to scaling in 
Eq. (|l9|) we have also included higher powers of q 3 L yi , which become particularly important when y\ is small (i.e. 
when a is close to d/2), higher powers of q\L d ~ 2y ^ = q\L~ d l 2 , and the crossterm proportional to L Vt+Vi . 

All analyses were carried out on the same data set as used in Ref. |5^, to which several data points have been added 
for most values of a. First, we have only kept fixed the exponents in the correction terms, t/i and y^. The corresponding 
estimates for Q and y* are shown in the third and fourth column of Table |. One observes that the Monte Carlo results 
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for both Q and y£ are in quite good agreement with the renormalization predictionsE3E3 Q = 87r 2 /r 4 (-|) = 0.456947. . . 
and yl — d/2. However, the uncertainties in the estimates increase considerably with increasing a, because the leading 
irrelevant exponent becomes very small. An exception is the relatively large uncertainty in y% {d = 1, a = 0.2), which 
originates from the fact that the Monte Carlo data were taken in a rather narrow temperature region around the 
critical point. Furthermore, an accurate simultaneous determination of Q and y* is very difficult, because of the 
correlation between the two quantities. Therefore we have repeated the same analysis with Q fixed at its theoretical 
prediction — as appears justified by the values for Q in Table | — in order to obtain more accurate estimates for y£. 
The results, shown in the fifth column of Table ||, are indeed in good agreement with the theoretically expected values 
(last column). Thus, we have kept the thermal exponent fixed at its theoretical value in the further analysis, just as 
in Ref. |53|. The corresponding results for Q and K c are shown in Table As discussed in Ref . |5^, over the full range 
of a and d the Monte Carlo results for Q show good agreement with the renormalization prediction, thus confirming 
the universality of this quantity above the upper critical dimension. In comparison with the estimates presented in 
Table I of Ref. [5^, two minor remarks apply. First, for Q(d = 3, a = 0.4) one decimal place too much was quoted, 
suggesting a too high accuracy. Secondly we note that the newest result for K c (d = 3, a = 1.2) deviates two standard 
deviations from the earlier estimate. 

The universality of Q is illustrated graphically in Figs. Q(a)-^(c), where the increasing importance of corrections 
to scaling upon approaching the upper critical dimension clearly follows from the size of the error bars. At the upper 
critical dimension itself (e = 0) this culminates in the appearance of logarithmic corrections, where the finite-size 
scaling form of Ql is given by 



Q L {T) = Q + viL yt {hxL) 



1/6 



t 



(InL) 2 / 3 



-paL^QnL) 1 / 3 



t + v- 



(lnL) 2 / 3 



+ qi L d - 2 '^ rr ■ ■■■ ■ (20) 

in Lj 

The ellipses denote terms containing higher powers of L d ~ 2yh and 1/lnL. The extremely slow convergence of this 
series is reflected in the uncertainty in the resulting estimates for Q at the upper critical dimension. To illustrate the 
dependence of the finite-size corrections on e more directly, Fig. ^(a) displays (for various values of a) the finite-size 
scaling functions as they follow from a least-squares fit of the data for d = 1 to Eqs. ( |l9| ) and (p0|), respectively. 
Although one clearly observes the increase of finite-size corrections when a — ► d/2, the true nature of the logarithmic 
corrections in ( p0| ) cannot be appreciated from this graph. To emphasize the difference between e — and e < 0, we 
therefore also show [Fig. §(b)] the same plot for the enormous range < L < 10 10 . Now it is evident how strongly 
the case e = differs even from a case with strong power-law corrections, such as a = 0.4 (e = —0.2). 

We have used the universality of Q to considerably narrow the error margins on K c by fixing Q at its theoretical 
value in the least-squares fit. The corresponding couplings are shown in Table |n| as well. The relative accuracy of the 
critical couplings lies between 1.5 x 10 -5 and 5.0 x 10~ 5 . For the one-dimensional case, we can compare these results 



to earlier estimates, see Table III. One notes that the_aewest estimates are more than two orders of magnitude more 
accurate than previous estimates. The first estimatesEa were obtained by carrying out exact calculations for chains 
of 1 to 20 spins and subsequently extrapolating these results using Pade approximants. Note that the estimates 
for T c in Ref. |l8| are expressed in units of the inverse of the Riemann zeta function and thus must be multiplied 
by £(l|=ker). All couplings are somewhat too high, but still in fair agreement with our estimates. The results of 
Domanlij have no error bars. Still, his results are worrying, since he carries out a cluster approach, obtaining critical 
couplings which start at the mean-field value for cluster size zero and increase monotonically with increasing cluster 
size, as they should, since mean- field theory yields a lower bound on the critical couplings (see below). Thus, he 
argues that the true couplings will lie higher than his best estimates (obtained for cluster size 10). However, all these 
best estimates lie already above our estimates, which seems to indicate a problem inherent in his approach. Ref. |2^ 
presents results of an approximation coined "finite-range scaling" with error margins of 1%. For a = 0.1 the error 
is considerably underestimated, but for the other values of the decay parameter the couplings agree with our results 
well within the quoted errors. The same technique was applied in Ref. but now the uncertainty in the couplings 
was estimated to be less than 10%, for small a a few times larger. This is clearly a too conservative estimate, as 



the difference with our results is only a few percent for a — 0.1 and considerably less for larger a. In Ref. 21, the 
coherent- anomaly method was used to obtain two different estimates without error margins. We have quoted the 
average of the two results, with their difference as a crude measure for the uncertainty. The agreement is quite good, 
although all results lie systematically above our values. Yet another approach has been formulated in Ref. where 
the Onsager reaction-field theory was applied to obtain a general expression for the critical coupling, 

r(l + <r)sin(7r<7/2) 

(1 — <7)7T- L + cr 
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Unfortunately, no estimate for the accuracy of this expression is given, but it seems to generally underestimate the 
critical coupling by a few percent- Finally, some estimates have recently been obtained by means of the real-space 
renormalization-group technique.!^ 



In addition, Monroe has calculated various bounds on the critical couplings as shown in Table IV. The Bethe lattice 
approximatiorJHJ was used to obtain both upper and lower bounds, to which our results indeed conform, although 
it must be said thaL-the upper bounds do not constitute a very stringent criterion. Furthermore, the application of 
Vigfusson's methodcJ has yielded even closer lower bounds for a — 0.1 and a — 0.2. 

Apart from these approximations, one may also use mean-sfield theory to make some predictions concerning the 
critical coupling in the limit a \ 0. It was shown by BrankovEj that in this limit the d-dimensional system with an 
interaction potential oc ajr^y" is equivalent to the Husimi-Temperley mean spherical model. More specifically, it 
was conjectured by CannascB that for the one-dimensional case lim CT ^o K c ~ ct/2, which is also the first term in the 
Taylor expansion of Eq. (pl|). Indeed, in mean- field theory one has zK^ F = 1, where z is the coordination number. 
For d — 1 this corresponds to the requirement 

2Kf F {a) £ ^ = 2iC IF (<x)C(l + = 1 , (22) 

71=1 

where denotes the Riemann zeta function. The expansion of C( x ) around x = 1 yields the conjectured relation 
lmvjo K^ F = c/2. Figure ^(a) shows the critical coupling as a function of the decay parameter a along with K^ AF (<r) 
and the asymptotic behavior for a J, 0. One observes that K c (a) indeed approaches K^ F (a) when a approaches 
zero. Furthermore, K^ F (a) is smaller than K c (a) for all a, as one expects from the fact that mean-field theory 
overestimates the critical temperature. It is interes ting to note that for a = 0.1 (K^ 0.047239) this lower bound 
already excludes the estimates given in Refs. ^ and ^] (cf. Table III). Replacing zK^ F by the integrated interaction, 
we can generalize such estimates to higher dimensionalities, 

n d/2 roo 1 

^ MF(a) rW dr ^ = 1 - (23) 

1 \2) Jm o 

For d > 1, the lower distance cutoff mo of the integral, i.e. the minimal interaction distance with the nearest neigh- 
bors, does not have an isotropic value, since there is no interaction within an elementary cube around the origin. 
Nevertheless, a constant value mo, e.g. mo = 1/2, is a good approximation. Furthermore, for d = 1 the integral is 
only a first-order approximation of Eq. (p2]), but for d = 2 and d — 3 it precisely corresponds to the interaction ( |lg| ) 
and its generalization to d = 3, respectively. As a first estimate one thus obtains 

\™Kf F {a) = V -%am° Q . (24) 

An expansion in terms of a shows that the first term is independent of mo. For d = 1,2,3 one finds, respectively, 
K™ F - cr/2, Kf F - o_A27r), K™ f - 0-/(4?!-). Figures|(b) and|(c) show K c (a) for d = 2 and d = 3, the corresponding 
asymptotes and Eq. (E4j) with mo = 1/2. 

The deviation of K c (<r) from K^ F (a) is also expressed by the last term in the renormalization expression @. 
However, in order to assess the a dependence of this term one has to calculate the a dependence of the coefficients a 
and c, arising from the integrals over the cr-dependent propagators. 

C. Determination of critical exponents 

1. Magnetic susceptibility 

The magnetic susceptibility x is directly proportional to the average square magnetization density, 

X = L d (m 2 ) , (25) 
and thus we can use Eq. (|lj) to analyze the finite-size data. Expanding this equation in t and u we obtain for e < 

X = L 2y *- d (ao + aifcZX + a 2 t 2 L 2 < +■■■ + hL n + • ■ ■) (26) 

and for e = 
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L 2v *- d VhTL 

x a + a 1 L yt (lnL) 1 / 6 



'(InL) 2 / 3 



+ a 2 L 2vt (In L) 1/3 [t + v 



(InL) 2 / 3 



InL 



(27) 



The analytic part of the free energy might give rise to an additional constant, but this could not be observed in our 
simulations, because it is dominated by the corrections to scaling. In Table ^ we list the results of an analysis of 
the numerical data. For all examined systems we have determined the exponent z/jj and the critical coupling. The 
estimates for the latter are in good agreement with those obtained from the analysis of the universal amplitude ratio Q. 
Furthermore, the exponents agree nicely, for all dimensionalities, with the renormalization prediction y^ — 3e?/4. Just 
as before, the uncertainties increase with increasing er, although the analyses at the upper critical dimension itself 
seem to yield better results than those just above it. Compare in particular the results for a = 1.4 (y; = —0.2) and 
(7 = 1.5. The logarithmic prefactor in Eq. (^) can be clearly observed in the sense that the quality of the least-squares 
fit decreases considerably when this factor is omitted. To reduce the uncertainty in the exponents we have repeated 
the analysis with K c fixed at the best values in Table [0|, i.e. those obtained with fixed Q. The corresponding estimates 
of yjj are also shown in Table [v| and are indeed in good agreement with the renormalization predictions. 

Now we can calculate the critical exponents and compare them to earlier estimates for d = 1. We do this for the 
correlation length exponent v = l/(y* +y;/2) and the magnetization exponent (3 — (d—y^)/y£. The results are shown 



in Tables VI and VII. Since all our estimates for yj* and j/jj agree with the renormalization values, also v and j3 are in 
agreement with the classical critical exponents. Unfortunately, the accuracy in both exponents is seriously hampered 
by the uncertainty in yl , which has only been determined from the temperature-dependent term in Q. In particular 
the results for v from Ref. [l2] are, for small er, in better agreement with the theoretically predicted values than our 
estimates. However, all previous results, both for v and for (3, deviate seriously from the predicted values when a 
approaches 1/2, which is not the case for our values. This can probably be attributed to the fact that corrections to 
scaling have been taken into account more adequately. 



2. Spin-spin correlation function 



In Sec. Ill two different decay modes for the spin-spin correlation function were derived. The relative magnitude of 
r and L determines which of the modes applies. In the bulk of our simulations we have restricted r in g(r) to r = L/2. 
Since this quantity reflects the k = mode of the correlation function, we write for e < an expression analogous to 
that for the magnetic susceptibility, 



g(L/2) = tfyt- 



2d 



c + citL Vi + c 2 



p L 2y' t 



d Y V>^ 



(28) 



and for e = 
g(L/2) = L 2vh - 2d VhTL 



c + 01^ (InL) 1 / 6 [t + v 



(InL) 2 / 3 



+ c 2 L 22/t (lnL) 1 / 3 [t + 



L~ Vt 
'(InL) 2 / 3 



di 
InL 



(29) 



For values of r such that g(r) does not correspond to this mode of the correlation function, the er-dependent exponent 
2/h w iU appear in (28) instead of y^. Furthermore, the logarithmic prefactor in ( p9|) will be ab sent, as it arises from 
the dangerous irrelevant variable [cf. Eq. (|Tq)]. The results of our analysis are shown in Table VIII. They evidently 
corroborate that the exponent y^ coincides with that appearing in the susceptibility. Also the factor VhiL in ( p9| ) 
was clearly visible in the least-squares analysis. The critical couplings agree with the estimates from Q and \ and we 
have again tried to increase the accuracy in y^ by repeating the analysis with K c fixed at their best values in Table ||. 
The accuracy of the results is somewhat less than of those obtained from the magnetic susceptibility, because we have 
now only used numerical data for even system sizes. The fact that the L dependence of g(L/2) is determined by the 
k = mode raises the question whether one can also observe the power-law decay described by ?/ in finite systems. To 
this end, we have sampled g(r) as a function of r in the one-dimensional model. In order to clearly distinguish the two 
predictions for the decay of g(r) we have examined a system far from the upper critical dimension, viz. with a = 0.1. 
It turned out to be necessary to sample very large system sizes to observe the regime where g(r) cx r _ ( d ~°'). Figure || 
displays the spin-spin correlation function scaled with L d / 2 versus r/L. The scaling makes the results collapse for 
r of the order of the system size. Here, the correlation function levels off. This is the mean-field like contribution 
to the correlation function, which dominates in the spatial integral yielding the magnetic susceptibility. For small r 
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the data do not collapse at all, which shows that g(r) exhibits different scaling behavior in this regime. Indeed, the 
correlation function decays here as r~( d ~ CT ) = 7-~ - 9 and not as r~ d / 2 . Note, however, that this regime is restricted to 
a small region of r and can only be observed for very large system sizes. 

It is interesting to note that already Nagle and BonneiuJ have tried to calculate rj in a spin chain with long- 
range interactions from finite-size data for the susceptibility. Because this calculation relied on the assumption that 
X(L, K c ) — \{L — 1, K c ) ~ g{L) ~ L _ ( d_2 +'w ! they called the corresponding exponent 77. The results for 77 turned out 
to assume a constant value approximately equal to 1.50 for < a < 0.5. Thus, the identification of fj with 77 was 
assumed to be invalid in Ref. |3f| Now we see that 77 is in excellent agreement with 77* = d + 2 — 2y£ — 2 — d/2. 

V. CONCLUSIONS 

In this paper we have studied systems with long-range interactions decaying as 7- _ ( d +°") m one, two, and three 
dimensions in the regime where these interactions exhibit classical critical behavior, i.e., for < a < d/2. From 
the renormalization equations we have derived the scaling behavior, including the corrections to scaling, for various 
quantities. These predictions, in particular the critical exponents and the scaling behavior of the amplitude ratio 
(m 2 ) 2 /(m 4 ), have been verified by accurate Monte Carlo results. At the upper critical dimension, the logarithmic 
factors appearing in the finite-size scaling functions could be accurately observed. The Monte Carlo results have been 
obtained with a dedicated algorithm. This algorithm is many orders of magnitude faster (up to the order of 10 s for 
the largest examined system) than a conventional Monte Carlo algorithm for these systems. Our analysis has also 
yielded estimates for the critical couplings. For d = 1 these values have an accuracy which is more than two orders of 
magnitude better than previous estimates and could thus serve as a check for half a dozen different approximations. 
For d — 2 and d — 3 we have, to our best knowledge, obtained the first estimates for the critical couplings. Finally, 
we have given both theoretical and numerical arguments that above the upper critical dimension the decay of the 
critical spin-spin correlation function in finite systems consists of two regimes: One where it decays as r _ ( c * _2 + ? ?) and 
one where it is independent of the distance. 

As an outlook we note that many interesting results may be expected below the upper critical dimension, where 
neither any rigorous-results nor any accurate numerical results are available. This regime will be the subject of a 
future investigation^ 

APPENDIX A: DETAILS OF THE MONTE CARLO ALGORITHM FOR LONG-RANGE 

INTERACTIONS 

The cluster algorithm applied in this study has been described for the first time in Ref. |34| A somewhat more 
elaborate treatment of the mathematical aspects was given in Ref. |6^. Although conceptually no new aspects arise in 
the application to algebraically decaying interactions in more than one dimension, several important practical issues 
must be taken care of in actual simulations. It is the purpose of this Appendix to discuss these issues and their 
solutions in some more detail. We do not repeat the full cluster algorithm here, but only describe how the cluster 
formation process proceeds from a given spin s, which has already been added to the cluster (the so-called current 
spin). 

As explained in Ref. [54|, the key element of the algorithm lies in splitting up the so-called bond-activation probability 
p(si,Sj) — 8 SiS pij — 5 SiS [l — exp(— 2Jij)] into two parts, namely the Kronecker delta testing whether the spins Si 
and Sj are parallel and the "provisional" bond-activation probability Pij . This enables us to define a cumulative bond 
probability C(fc), from which we can read off which bond is the next one to be provisionally activated, 

C(j) = £p(n) (Al) 

n=l 

with 

P(Tl) = (1 - P l)(l - P2 ) ■ ■ ■ (1 - Pn-l)Pn • (A2) 

Pj = 1 — exp(— 2Jj) is an abbreviation for ppj, i.e., we define the origin at the position of the current spin. When 
comparing the expressions to those in Ref. |5J one must take into account that we now are working with Ising instead of 
Potts couplings. P{n) is the probability that in the first step 77 — 1 bonds arc skipped and the nth bond is provisionally 
activated. Now the next bond j that is provisionally activated is determined by a random number g £ [0, 1): j — 1 
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bonds are skipped if C(j — 1) < <7 < C(j). The number j can be easily determined from g once we have tabulated the 
quantity C(j) in a look-up table. If the jth bond is placed to a spin Sj that is indeed parallel to the current spin s, 
then Sj is added to the cluster (i.e., the jth bond is activated). Subsequently we skip again a number of bonds before 
another bond at a dista nce k > j is provisionally activated. The appropriate cumulative probability is now given by 
a generalization of Eq. (|AlD (see Ref. p4) , 



n=j+l 



H 0--Pm) 



m—j-\-l 



p n = l- exp -2 ^ Jn 



(A3) 



n=j+l 



In principle we need now for each value of j another look-up ta ble contai ning the Cj(k). This is hardly feasible and 
fortunately not necessary, as follows from a comparison of Eqs. (Al) and (A3). Namely, 



C(k) = C (k) = C{j) + 



11(1 -Pi) 



C j (k) = C(j) + [l-C(j)]C j (k) 



(A4) 



or Cj(k) = [C(k) — C(j)]/[1 — C(j)]. So we can calculate Cj(k) directly from C(k). In practice one realizes this by 
using the bond distance j of the previous bond that was provisionally activated to rescale the (new) random number 
g to g' S [C(j),l); g' — C(j) + [1 — C(j)]g. Since we consider only ferromagnetic interactions, lim^oo C(j) exists 
and is smaller than 1, cf. Eq. (A3). Still we can accommodate only a limited number of bond distances in our look-up 
table and must therefore devise some approximation scheme to handle the tail of the long-range interaction, which is 
essential for the critical behavior. This issue is addressed below. Furthermore, this description only takes into account 
the bonds placed in one direction. The actual implementation of the algorithm must of course allow for bonds in both 
directions (assuming that d = 1). 

An alternative for the look-up table exists for interactions which can be explicitly summed. In those cases, Eq. ( |A3| ) 
can be solved for k, yielding an expression for the bond distance in terms of Cj(fc), i.e., in ter ms of the random 
number g. For the interaction defined in Sec. || the sum appearing in the right-hand side of (A3) is (for j = 0) the 
truncated Riemann zeta function, 



,d+ty 



(A5) 



71=1 



which cannot be expressed in closed form. In more than one dimension, a look-up table is very impractical and an 
interaction which can be summed explicitly becomes very desirable. Therefore we have taken an isotropic, continuous 
interaction of the form J = K/r d+a . The interaction with a spin at lattice site n is then given by the integral of J 
over the elementary square (cube) centered around n [cf. Eq. (pf)] and the cumulative bond probability yields the 
(not necessarily integer-valued) distance k at which the first provisional bond is placed. To this end, the sum in ( |A3| ) 
is replaced by a ci-dimensional integral over the coupling J . As J is isotrop ic, o nly an integral over the radius remains, 
which runs from the minimal bond distance up to k. Thus for d = 2 Eq. (A3) reduces to 



Cj(k) = 1 — exp 



2<kK ( 1 



1 



3 



and in d = 3 the factor 27r is simply replaced by Ait. Equating Cj(k) to the random number g we find 



-l/a 



(A6) 



(A7) 



Rescaling of the random number is no longer required: The lowest value, g = 0, leads to a provisional bond at the 
same distance as the previous one, k — j. If g = Cj(oo) = 1 — exp[— (AttK / o-)j~ a ] the next provisional bond lies at 
infinity and thus g G [Cj(oo), 1) yields no bond at all. Once the distance k has been obtained, d — 1 further random 
numbers g\, . . . are required to determine the direction of the bond. In d = 2, we set 4> = g\/ (27r). The coordinates 
of the next provisional bond (relative to the current spin) are then (r x ,r y ) = (k cos <f>, k sin 0), which are rounded 
to the nearest integer coordinates. Finally, the periodic boundary conditions are applied to map these coordinates 
onto a lattice site. For the next provisional bond, j is set equal to k {not to the rounded distance!) and a new k is 
determined. If no bond has been placed yet, j is set to 1/2, the lowest possible bond distance. Hence it is possible 
to find a 1/2 < k < V2/2 and an angle 4> such that the corresponding lattice site is the origin. This does not affect 
the bond probabilities, but it is of course a "wasted" Monte Carlo step. For d = 3 the process is similar, except that 



11 



we need another random number gi to determine a second angle — tt/2 < tp < 7r/2, such that sin?/' is distributed 
uniformly; sin^ = 1 — 2g2- The bond coordinates are given by (A; cos ?/> cos (/>, kcosipsmcj), ksiatp). 

This approach can also be applied in the one-dimensional case, where the geometrical factor 2ir in ( |A6| ) must be 
replaced by 2, which reflects the fact that bonds can be put to the left and to the right of the origin. The direction 
of the bond is then simply determined by another random number. As has already been mentioned in Ref. p4l t his 
can be used to cope with the limited size M of the look-up table. Beyond the bond distance M the sum in ( |A3| ) is 
approximated by an integral. I.e., if the random num ber g lies in the interval [C(M), C(oo)), the bond distance k is 
determined from the one-dimensional version of (|A7| ) , where the lower part of the integral is replaced by an explicit 
sum 



Here, the geometrical factor is absent, as we have opted to treat "left" and "right" separately in our simulations (no 
additional random number is required in that case). The approximation ( |A8| ) effectively introduces a modification of 
the spin-spin interaction, which however can be made arbitrarily small by increasing M. Note that the offset 1/2 in 
the first term ensures a precise matching of the discrete sum and the integral approximation: the random number 
g = C(M) = 1 — exp[— 2K^2^f = i n^^ 1+<J ^\ yields k = M + 1/2 which is precisely the lowest k that is rounded to the 
integer bond distance M + 1. 

The accuracy of this procedure is further limited by the finite resolution of random numbers. E.g., in our simulations 
the original random numbers are integers in the range [0, 2 32 - 1]. Thus, for bond distances I such that C(7) - C(l - 1) 
is of the order 2~ 32 , the discreteness of the random numbers is no longer negligible. For d = 2 and d — 3, the 
discreteness of the angles also limits the lattice sites that can be selected for a provisional bond, but this generally 
occurs at distances larger than I. Once the value of I has been determined, with a safe margin, there are various 
approaches to this limitation. One may, e.g., draw another random number to determine the precise bond distance. A 
simpler approach is to distribute all bonds beyond / uniformly over the lattice, in order to prevent that certain lattice 
sites are never selected. However, one should take care that such simple approaches do not essentially modify the 
critical behavior. If I is relatively small, the error introduced by a random distribution of the bond distances might be 
larger than the effect of an interaction which decreases slightly nonmonotonically at large distances. Furthermore, in 
order to preserve the symmetry of the lattice, such a uniform distribution of the bonds should occur outside a square 
(cube) instead of a circle (sphere) with radius I. 
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FIG. 1. The amplitude ratio Q as a function of the decay parameter a in (a) d = 1, (b) d = 2, and (c) d = 3 dimensions. 
The solid line marks the renormalization prediction. 
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FIG. 2. The amplitude ratio Q in a one- dimensional system as a function of the system size L for various values of a. 
Figure (a) illustrates the increase of the finite-size corrections when the upper critical dimension (a = d/2) is approached. 
Figure (b) emphasizes the difference between finite-size corrections above the upper critical dimension (power-law) and at the 
upper critical dimension itself (logarithmic). 
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FIG. 3. The critical coupling K c as a function of the decay parameter a for (a) d — 1, (b) d — 2, and (c) d = 3. Also 
shown is the asymptotic behavior for a J. as predicted by mean-field theory and mean-field values for K c over the full range 
of < a < d/2 (for d = 2 and d = 3 only approximately). 
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FIG. 4. The spin-spin correlation function versus r/L in the one-dimensional model with a = 0.1. Results for various system 
sizes are shown. For a discussion see the text. 
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TABLE I. The amplitude ratio Q and the thermal exponent j/t for systems with long-range interactions in one, two, and 
three dimensions, for several values of the decay parameter < a < d/2. The values in the fifth column have been obtained 
with Q fixed at the theoretically predicted value (see text) and the last column lists the renormalization predictions for yX . 



d 


a 


Q 


V* 


Vt 


RG 


1 


0.1 


0.4566 (8) 


0.507 (7) 


0.507 (7) 


i 

2 


1 


0.2 


0.455 (4) 


0.54 (4) 


0.504 (12) 


\ 
2 


1 


0.25 


0.457 (3) 


0.500 (8) 


0.500 (5) 


\ 

2 


1 


0.3 


0.454 (2) 


0.519 (14) 


0.506 (12) 


\ 

2 


1 


0.4 


0.457 (3) 


0.50 (2) 


0.50 (2) 


2 


1 


0.5 


0.462 (6) 


0.51 (5) 


0.49 (2) 


I 
2 


2 


0.2 


0.4574 (10) 


1.01 (2) 


1.01 (2) 


1 


2 


0.4 


0.455 (2) 


1.02 (2) 


1.009 (15) 


1 


2 


0.6 


0.450 (6) 


1.04 (4) 


1.008 (17) 


1 


2 


0.8 


0.454 (6) 


1.03 (9) 


1.03 (3) 


1 


2 


1.0 


0.450 (10) 


1.02 (3) 


1.03 (2) 


1 


3 


0.2 


0.4581 (11) 


1.51 (3) 


1.513 (18) 


2 


3 


0.4 


0.4561 (10) 


1.521 (18) 


1.512 (15) 


3 
2 


3 


0.6 


0.453 (3) 


1.53 (4) 


1.521 (14) 


3 
2 


3 


0.8 


0.458 (2) 


1.48 (2) 


1.487 (10) 


3 
2 


3 


1.0 


0.453 (10) 


1.52 (7) 


1.508 (9) 


3 
2 


3 


1.2 


0.447 (8) 


1.56 (2) 


1.519 (10) 


3 
2 


3 


1.4 


0.454 (5) 


1.48 (3) 


1.48 (3) 


3 
2 


3 


1.5 


0.449 (8) 


1.53 (5) 


1.46 (3) 


3 
2 



TABLE II. The amplitude ratio Q and critical couplings K c for systems with long-range interactions in one, two, and three 
dimensions, for several values of the decay parameter < a < d/2. The thermal exponent (see Table |l|) was kept fixed at its 
theoretical value in all analyses. The estimates for K c in the last column have been obtained by fixing Q at its renormalization 
prediction. The numbers between parentheses represent the errors in the last decimal places. 



d 


a 


Q 


K c 


K c 


1 


0.1 


0.4565 (8) 


0.0476162 (13) 


0.0476168 (6) 


1 


0.2 


0.4579 (7) 


0.092234 (2) 


0.0922314 (15) 


1 


0.25 


0.4579 (15) 


0.114143 (4) 


0.1141417 (19) 


1 


0.3 


0.4567 (15) 


0.136113 (4) 


0.136110 (2) 


1 


0.4 


0.457 (3) 


0.181151 (8) 


0.181150 (3) 


1 


0.5 


0.463 (5) 


0.229157 (8) 


0.229155 (6) 


2 


0.2 


0.4573 (10) 


0.028533 (3) 


0.0285324 (14) 


2 


0.4 


0.4565 (17) 


0.051824 (4) 


0.0518249 (14) 


2 


0.6 


0.456 (4) 


0.071364 (7) 


0.071366 (2) 


2 


0.8 


0.458 (5) 


0.088094 (7) 


0.088094 (2) 


2 


1.0 


0.447 (8) 


0.102556 (5) 


0.102558 (5) 


3 


0.2 


0.4584 (9) 


0.0144361 (10) 


0.0144354 (6) 


3 


0.4 


0.4569 (8) 


0.0262927 (16) 


0.0262929 (7) 


3 


0.6 


0.4581 (9) 


0.036050 (2) 


0.0360469 (11) 


3 


0.8 


0.4562 (13) 


0.044034 (2) 


0.0440354 (10) 


3 


1.0 


0.4571 (14) 


0.050515 (2) 


0.0505152 (12) 


3 


1.2 


0.457 (3) 


0.055682 (3) 


0.0556825 (14) 


3 


1.4 


0.455 (5) 


0.059666 (2) 


0.0596669 (11) 


3 


1.5 


0.449 (7) 


0.061251 (2) 


0.061253 (2) 



18 



TABLE III. Comparison between our best estimates of the critical couplings K c for the one-dimensional system and earlier 
estimates. 



a 


This work 


Ref. li 


Ref. |T| 


Ref. pc] 


Ref. P 


Ref. pi 


Ref. p7| 


Ref. p| 


0.1 


0.0476168 (6) 




0.0478468 


0.0505 (5) 


0.04635 


0.04777 (12) 


0.0469 


0.0481 


0.2 


0.0922314 (15) 


0.0926 (5) 


0.0933992 


0.0923 (9) 


0.09155 


0.0928 (3) 


0.0898 




0.25 


0.1141417 (19) 












0.1106 




0.3 


0.136110 (2) 


0.1370 (7) 


0.138478 


0.1362 (14) 


0.1359 


0.1375 (10) 


0.1314 


0.144 


0.4 


0.181150 (3) 


0.1825 (10) 


0.184081 


0.1815 (18) 


0.1813 


0.183 (2) 


0.1750 




0.5 


0.229155 (6) 


0.2307 (14) 


0.230821 


0.230 (2) 


0.2295 


0.231 (4) 


0.2251 


0.250 



TABLE IV. Comparison of our best estimates of the critical couplings for the one-dimensional system with some lower and 
upper bounds. 



a 


This work 


Ref. 


24 




Ref. 24 




Ref. 


25 




0.1 
0.2 
0.3 
0.4 
0.5 


0.0476168 (6) 
0.0922314 (15) 
0.136110 (2) 
0.181150 (3) 
0.229155 (6) 


> 0.04726 

> 0.08947 

> 0.1273 

> 0.1615 

> 0.1923 


< 0.09456 

< 0.1792 

< 0.2558 

< 0.3258 

< 0.3903 


> 0.04753 

> 0.09162 



TABLE V. Estimates for the critical coupling K c and the exponent y^ as obtained from the analysis of the magnetic 



susceptibility. The values for y^ in the fifth column have been obtained by fixing K c at their best estimates from Table |lf the 
error margins do not include the uncertainty in these values for K c . 



d 


a 


K c 


Vh 


Vh 


RG 


1 


0.1 


0.0476161 (19) 


0.7487 (14) 


0.7493 (6) 


3 

4 


1 


0.2 


0.092239 (4) 


0.752 (2) 


0.7504 (10) 


3 

4 


1 


0.25 


0.114145 (4) 


0.7477 (15) 


0.7468 (16) 


3 

4 


1 


0.3 


0.136110 (5) 


0.747 (3) 


0.7490 (17) 


3 

4 


1 


0.4 


0.181170 (10) 


0.749 (5) 


0.746 (3) 


3 

4 


1 


0.5 


0.229153 (6) 


0.748 (2) 


0.7490 (8) 


3 
4 


2 


0.2 


0.028537 (5) 


1.500 (6) 


1.495 (3) 


3 
2 


2 


0.4 


0.051830 (6) 


1.498 (9) 


1.496 (3) 


3 
2 


2 


0.6 


0.071370 (5) 


1.497 (6) 


1.498 (2) 


3 
2 


2 


0.8 


0.088095 (10) 


1.496 (5) 


1.495 (3) 


3 
2 


2 


1.0 


0.102556 (3) 


1.495 (4) 


1.497 (2) 


3 
2 


3 


0.2 


0.0144347 (9) 


2.249 (2) 


2.2504 (8) 


9 
4 


3 


0.4 


0.026296 (2) 


2.250 (6) 


2.246 (3) 


9 

4 


3 


0.6 


0.036046 (3) 


2.246 (7) 


2.244 (5) 


9 

4 


3 


0.8 


0.0440349 (17) 


2.243 (4) 


2.246 (3) 


9 

4 


3 


1.0 


0.050516 (3) 


2.239 (2) 


2.243 (7) 


9 

4 


3 


1.2 


0.055679 (2) 


2.247 (11) 


2.251 (7) 


9 
4 


3 


1.4 


0.0596636 (18) 


2.27 (3) 


2.26 (2) 


9 

4 


3 


1.5 


0.061251 (2) 


2.257 (12) 


2.249 (5) 


9 
4 
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TABLE VI. The correlation length exponent v as a function of a for the one-dimensional model, together with earlier 
estimates and the renormalization predictions. 





TniQ "wriTk" 

_1_ lllO WU1J\ 


Ref. 20 


Ref. 42 


Ref. 43 


RG 


0.1 


9.3 (6) 


9.12 


9.9 


10.48 


10.0 


0.2 


4.9 (3) 


4.90 


4.95 




5.0 


0.25 


4.00 (8) 








4.0 


0.3 


3.27 (12) 


3.41 


3.32 


3.90 


3.3... 


0.4 


2.50 (13) 


2.71 


2.68 




2.5 


0.5 


2.04 (8) 


2.34 


2.33 


2.81 


2.0 



TABLE VII. The magnetization exponent (5 as a function of a for the one-dimensional model, together with earlier estimates 
and the renormalization predictions. 



a 


This work 


Ref. |l| 


Ref. |l] 


RG 


0.1 

0.2 

0.25 

0.3 

0.4 

0.5 


0.494 (8) 
0.495 (13) 
0.506 (8) 
0.497 (15) 
0.51 (2) 
0.51 (2) 


0.5 

0.48 
0.45 
0.39 


0.495 
0.482 

0.460 
0.435 
0.408 


i 

2 
1 
2 
1 
2 
1 
2 
1 
2 
1 
2 



TABLE VIII. Estimates for the critical coupling K c and the exponent y^ as obtained from the analysis of the spin-spin 



correlation function. The values for y^ in the fifth column have been obtained by fixing K c at their best estimates from Table |l| 
the error margins do not include the uncertainty in these values for K c . 



d 


a 


K c 


Vh 


2/h 


RG 


1 


0.1 


0.047619 (3) 


0.750 (2) 


0.7488 (9) 


3 

4 


1 


0.2 


0.092233 (7) 


0.749 (3) 


0.7513 (16) 


3 

4 


1 


0.25 


0.114148 (10) 


0.750 (5) 


0.747 (2) 


3 

4 


1 


0.3 


0.136116 (7) 


0.753 (5) 


0.752 (3) 


3 

4 


1 


0.4 


0.181158 (15) 


0.747 (7) 


0.750 (4) 


3 

4 


1 


0.5 


0.229150 (7) 


0.749 (2) 


0.7503 (10) 


3 

4 


2 


0.2 


0.028535 (7) 


1.499 (9) 


1.496 (3) 


3 
2 


2 


0.4 


0.051831 (6) 


1.505 (6) 


1.499 (4) 


3 
2 


2 


0.6 


0.071369 (6) 


1.507 (4) 


1.502 (4) 


3 
2 


2 


0.8 


0.088091 (6) 


1.495 (7) 


1.497 (3) 


3 
2 


2 


1.0 


0.102554 (4) 


1.490 (6) 


1.496 (3) 


3 
2 


3 


0.2 


0.0144348 (16) 


2.256 (6) 


2.254 (4) 


9 

4 


3 


0.4 


0.026296 (3) 


2.257 (8) 


2.245 (5) 


9 

4 


3 


0.6 


0.036053 (4) 


2.262 (10) 


2.246 (4) 


9 

4 


3 


0.8 


0.044035 (4) 


2.252 (11) 


2.250 (5) 


9 

4 


3 


1.0 


0.050511 (5) 


2.228 (15) 


2.249 (9) 


9 

4 


3 


1.2 


0.055680 (3) 


2.253 (14) 


2.257 (9) 


9 

4 


3 


1.4 


0.059667 (2) 


2.22 (4) 


2.31 (4) 


9 

4 


3 


1.5 


0.061251 (5) 


2.26 (3) 


2.248 (7) 


9 

4 
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